5  Diagnósticos de Residuos

5.1 Pretratamientos de la Información

Cargar Recursos Externos

source("scripts/recursos.R")
source("scripts/funciones.R")

Definición de variables generales

## geographic projections
crs_latlon <- "+proj=longlat +datum=WGS84 +no_defs"
crs_utm <- 
  "+proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

Definición de variables específicas

# List of crimes/delitos and cities/ciudades
crimes = c(
  "Minor property offenses",
  "Robberies",
  "Burglaries",
  "Drunkennes, damages and disorders",
  "Family violence and aggression",
  "Injuries, drugs and weapons",
  "High violence and murder"
)
delitos = c("hurto", "robo", "asalto", "amen", "abus", "crim", "viol")

#ciudad de ejemplo
cities <- c("Coquimbo-Serena")
ciudades <- c("urb_4_18")
obs <- 1
delito <- delitos[1]


Lectura de Datos Espaciales

Para el caso del presente ejemplo práctico se utiizará la ciudad de Coquimbo-Serena.

city <- ciudades[obs]
cityp <- readRDS(paste0("data/",city,".rds"))
cityp <- spTransform(cityp, CRS(crs_utm))

El objeto espacial cityp cuyos puntos corresponde donde ha ucurrido eventos deictuales, cuyo formato SpatialPointDataFrame.


Matriz de vecindad espacial: Se utilizó la función propia spatial_weights() @(sec-spatial-weights)

# matriz de vecindad espacial
nb <- spatial_weights(city_df = cityp, nvec= 12)


Transformaciones a los datos.

Se utilizó funciónes de transformacion general mencionadas en Section 6.2 que corresponden a las siguientes:

  • Transformación Cúbica
  • Transformación Logarítmica
  • Cáculo del Lag Espacial
cityp <- cityp %>%
    st_as_sf(crs = 4326) %>% 
    mutate(across(.cols = all_of(delitos),
                  .fns = t_cub, .names = "cub{.col}")) %>% 
    mutate(across(.cols = all_of(delitos),
                  .fns = t_log1p, .names = "log{.col}")) %>% 
    mutate(across(.cols = all_of(delitos),
                  .fns = lag_spatial, nb, .names = "lag{.col}")) %>%
    as("Spatial")
zona id x y pflot viol crim abus amen robo hurto asalto denspob educ tamhog constr ofcom vulzon mobil lpflot neduc ldenspob lofcom lconstr cubhurto cubrobo cubasalto cubamen cubabus cubcrim cubviol loghurto logrobo logasalto logamen logabus logcrim logviol laghurto lagrobo lagasalto lagamen lagabus lagcrim lagviol
4101161005 22 -71.239 -29.862 1.000 0 0 0 0 0 0 0 144.756 10.867 3.170 4328.618 0.000 0.367 0.079 0.178 0.477 0.767 0.000 0.840 0 0 0 0 0 0 0 0 0 0.000 0.000 0.000 0 0 0.083 0.083 0.083 0.417 0.250 0.083 0
4101161005 32 -71.243 -29.863 1.000 0 0 0 0 0 0 0 227.772 10.716 3.393 3830.009 0.000 0.367 0.079 0.178 0.466 0.846 0.000 0.829 0 0 0 0 0 0 0 0 0 0.000 0.000 0.000 0 0 0.000 0.083 0.083 0.917 0.583 0.083 0
4101161005 35 -71.240 -29.863 1.034 0 0 0 0 0 0 1 171.133 11.801 3.172 4328.618 0.000 0.367 0.079 0.182 0.550 0.796 0.000 0.840 0 0 1 0 0 0 0 0 0 0.693 0.000 0.000 0 0 0.000 0.083 0.000 0.500 0.250 0.083 0
4101161005 36 -71.239 -29.863 1.000 0 0 0 0 0 0 0 154.603 10.870 3.135 4328.618 0.000 0.367 0.079 0.178 0.478 0.779 0.000 0.840 0 0 0 0 0 0 0 0 0 0.000 0.000 0.000 0 0 0.083 0.083 0.083 0.417 0.250 0.083 0
4101161005 37 -71.238 -29.863 1.025 0 0 0 0 0 0 0 155.145 10.858 3.148 3230.309 0.000 0.367 0.079 0.181 0.477 0.779 0.000 0.815 0 0 0 0 0 0 0 0 0 0.000 0.000 0.000 0 0 0.083 0.083 0.083 0.000 0.250 0.083 0
4101161001 48 -71.244 -29.864 1.125 0 0 0 0 0 0 0 168.179 11.958 3.336 2384.222 14.682 0.292 0.122 0.193 0.562 0.793 0.574 0.788 0 0 0 0 0 0 0 0 0 0.000 0.000 0.000 0 0 0.000 0.000 0.250 0.500 0.417 0.083 0
4101161005 49 -71.243 -29.864 1.071 0 0 0 0 0 0 0 212.743 10.978 3.415 3795.759 0.000 0.367 0.079 0.187 0.486 0.834 0.000 0.828 0 0 0 0 0 0 0 0 0 0.000 0.000 0.000 0 0 0.000 0.083 0.250 0.917 0.583 0.083 0
4101161005 50 -71.242 -29.864 1.000 0 0 1 1 0 0 0 220.188 10.705 3.409 4026.613 0.000 0.367 0.079 0.178 0.465 0.840 0.000 0.834 0 0 0 1 1 0 0 0 0 0.000 0.693 0.693 0 0 0.000 0.083 0.000 1.000 0.417 0.083 0
4101161005 51 -71.241 -29.864 1.021 0 0 1 1 0 0 0 243.501 10.120 3.337 4328.618 0.000 0.367 0.079 0.180 0.419 0.857 0.000 0.840 0 0 0 1 1 0 0 0 0 0.000 0.693 0.693 0 0 0.000 0.083 0.083 0.667 0.250 0.000 0
4101161005 52 -71.240 -29.864 1.022 0 0 0 0 0 0 0 320.023 10.720 3.324 4328.618 0.000 0.367 0.079 0.180 0.466 0.905 0.000 0.840 0 0 0 0 0 0 0 0 0 0.000 0.000 0.000 0 0 0.000 0.083 0.167 0.583 0.250 0.000 0


5.2 Evaluación De Residuos

5.2.1 Modelos MLB

Lectura del Modelo.

Warning in spTransform(xSP, CRSobj, ...): NULL source CRS comment, falling back
to PROJ string
Warning in wkt(obj): CRS object has no comment
model_mlb <- readRDS(paste0("output/mlb_",city,"_",delito,".rds"))

5.2.1.1 Extracción de Residuos

cityp$res_fit_mlb <- residuals(model_mlb)
cityp$fitted_fit_mlb <- fitted(model_mlb)
cityp$sd_breaks_mlb <- scale(cityp$res_fit_mlb)[,1]
breaks_qt <- classInt::classIntervals(cityp$sd_breaks_mlb, 
                                      n = 5, style = "fisher")
my_breaks <- round(breaks_qt$brks,2)
city_sf <- cityp %>% st_as_sf(crs = 32719) 
residuos_mv <- mapview(city_sf,zcol="res_fit_mlb", at =my_breaks, cex=2) 
residuos_mv

5.2.1.2 Global Moran

Global Moran a Variable Depediente

spdep::moran.test(x = cityp$laghurto,  listw = nb)

    Moran I test under randomisation

data:  cityp$laghurto  
weights: nb    

Moran I statistic standard deviate = 174.4, p-value <
0.00000000000000022
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     0.7640052224     -0.0001215805      0.0000191963 
# moran.plot(x =cityp$laghurto,  listw = nb, labels=as.character(cityp$id))

Global Moran de Residuos

spdep::moran.test(x = cityp$res_fit_mlb,  listw = nb)

    Moran I test under randomisation

data:  cityp$res_fit_mlb  
weights: nb    

Moran I statistic standard deviate = 19.256, p-value <
0.00000000000000022
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
    0.08478830236    -0.00012158055     0.00001944476 
# moran.plot(x =cityp$res_fit_mlb,  listw = nb, labels=as.character(cityp$id))

5.2.1.3 Local Moran de Residuos

# Calcular Local Moran
lmoran = localmoran(cityp$res_fit_mlb, nb)
cityp$res_scaled = as.numeric(scale(cityp$res_fit_mlb))
cityp$lag_scaled = lag.listw(nb, cityp$res_scaled)
cityp@data = cbind(cityp@data, lmoran = as.data.frame(lmoran)[, 5])

# Umbral de significancia estadistica
pval=0.05

# Definir cuadrantes
cityp[(cityp$res_scaled >= 0 & cityp$lag_scaled >= 0) & (cityp$lmoran <= pval), "clusterM"] = "HH" # plot
cityp[(cityp$res_scaled <= 0 & cityp$lag_scaled <= 0) & (cityp$lmoran <= pval), "clusterM"] = "LL" # plot
cityp[(cityp$res_scaled >= 0 & cityp$lag_scaled <= 0) & (cityp$lmoran <= pval), "clusterM"] = "HL"
cityp[(cityp$res_scaled <= 0 & cityp$lag_scaled >= 0) & (cityp$lmoran <= pval), "clusterM"] = "LH"
cityp[(lmoran[, 5] > pval), "clusterM"] = "NS"
table(cityp$clusterM)

  HH   HL   LH   LL   NS 
 284   33  295  331 7283 
city_sf <- cityp %>% st_as_sf(crs = 32719) %>% 
  filter(clusterM != "NS")

residuos_moran <- mapview(city_sf,zcol="clusterM", cex=2 )
residuos_moran

5.2.2 Modelos SAR

Warning in spTransform(xSP, CRSobj, ...): NULL source CRS comment, falling back
to PROJ string
Warning in wkt(obj): CRS object has no comment
model_sar <- readRDS(paste0("output/sar_",city,"_",delito,".rds"))

5.2.2.1 Extracción de Residuos SAR

cityp$res_fit_sar <- residuals(model_sar)
cityp$fitted_fit_sar <- fitted(model_sar)
This method assumes the response is known - see manual page
cityp$sd_breaks_sar <- scale(cityp$res_fit_sar)[,1]
breaks_qt <- classInt::classIntervals(cityp$sd_breaks_sar, 
                                      n = 5, style = "fisher")
my_breaks <- round(breaks_qt$brks,2)
city_sf <- cityp %>% st_as_sf(crs = 32719) 

residuos_mv <- mapview(city_sf,zcol="res_fit_sar", at =my_breaks, cex=2) 
  
residuos_mv

5.2.2.2 Global Moran

Global Moran a Variable Depediente

spdep::moran.test(x = cityp$cubhurto,  listw = nb)

    Moran I test under randomisation

data:  cityp$cubhurto  
weights: nb    

Moran I statistic standard deviate = 79.737, p-value <
0.00000000000000022
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
    0.35153153375    -0.00012158055     0.00001944977 
# moran.plot(x =cityp$laghurto,  listw = nb, labels=as.character(cityp$id))

Global Moran de Residuos

spdep::moran.test(x = cityp$res_fit_sar,  listw = nb)

    Moran I test under randomisation

data:  cityp$res_fit_sar  
weights: nb    

Moran I statistic standard deviate = -1.2062, p-value = 0.8861
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
   -0.00544188551    -0.00012158055     0.00001945508 
# moran.plot(x =cityp$res_fit_mlb,  listw = nb, labels=as.character(cityp$id))

5.2.2.3 Local Moran de Residuos

# Calcular Local Moran
lmoran = localmoran(cityp$res_fit_sar, nb)

cityp$res_scaled = as.numeric(scale(cityp$res_fit_sar))
cityp$lag_scaled = lag.listw(nb, cityp$res_scaled)
cityp@data = cbind(cityp@data, lmoran = as.data.frame(lmoran)[, 5])

# Umbral de significancia estadistica
pval=0.05

# Definir cuadrantes
cityp[(cityp$res_scaled >= 0 & cityp$lag_scaled >= 0) & (cityp$lmoran <= pval), "clusterM"] = "HH" # plot
cityp[(cityp$res_scaled <= 0 & cityp$lag_scaled <= 0) & (cityp$lmoran <= pval), "clusterM"] = "LL" # plot
cityp[(cityp$res_scaled >= 0 & cityp$lag_scaled <= 0) & (cityp$lmoran <= pval), "clusterM"] = "HL"
cityp[(cityp$res_scaled <= 0 & cityp$lag_scaled >= 0) & (cityp$lmoran <= pval), "clusterM"] = "LH"
cityp[(lmoran[, 5] > pval), "clusterM"] = "NS"
table(cityp$clusterM)

  HH   HL   LH   LL   NS 
 214   19  156   32 7805 
city_sf <- cityp %>% st_as_sf(crs = 32719) %>% 
  filter(clusterM != "NS")

residuos_moran <- mapview(city_sf,zcol="clusterM", cex=2 )
residuos_moran